May 8, 2017

Problema modelo

Se mide una variable en una población que viene de dos subpoblaciones normales independientes, esta variable se distribuye en cada subpoblación con media y varianza diferentes. El resultado de esto es una mixtura.

El objetivo es estimar la media y la varianza de las subpoblaciónes, asi como la probabilidad de que un individuo pertenezca a alguna de las dos poblaciones.

Es importante notar que no se tiene información de a que subpoblación pertenece cada individuo.

Población

Módelo

La densidad del módelo es:

\[ f(x;\mu_1,\mu_2,\sigma^2_1,\sigma^2_2,\pi)=\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2) \]

Con esto la función de log-verosimilitud queda:

\[ l(x;\mu_1,\mu_2,\sigma^2_1,\sigma^2_2,\pi)= \]

\[\sum_{i=1}^n \ln(\pi d(x_i;\mu_1,\sigma^2_1)+(1-\pi)d(x_i;\mu_2,\sigma^2_2)) \]

Ecuaciónes normales

Obteniendo un sistema a partir de las derivadas tenemos:

\[ \frac{dl}{d\pi}=\sum_{i=1}^n \frac{d(x_i;\mu_1,\sigma^2_1)-d(x_i;\mu_2,\sigma^2_2)}{\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2)} \] \[ \frac{dl}{d\mu_k}=\sum_{i=1}^n \frac{(-1)^{k+1}((3-2k)\pi+k-1)\frac{d}{d\mu_k}d(x_i;\mu_k,\sigma^2_k)}{\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2)} \]

\[ \frac{dl}{d\sigma^2_k}=\sum_{i=1}^n \frac{(-1)^{k+1}((3-2k)\pi+k-1) \frac{d}{d\sigma^2_k}d(x_i;\mu_k,\sigma^2_k)}{\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2)} \]

El método de Newton

Sea \(f(x)\) una función de la cúal deseamos encontrar una raiz, sea \(\alpha\) esta raiz. Si realizamos la expansión de Taylor en el punto \(x_i\)

\[ f(x)=f(x_i)+f'(x_i)(x-x_i)+O((x-x_i)^2) \]

Si evaluamos en \(\alpha\) e ignoramos el segundo termino,

\[ 0=f(\alpha)=f(x_i)+f'(x_i)(\alpha-x_i) \] Así \[ \alpha=x_i-\frac{f(x_i)}{f'(x_i)} \]

Método de Newton

Este \(\alpha\) que estamos proponiendo no es exactamente la raiz (ya que aproximamos al descartar los terminos de orden 2), pero esperamos que este más cerca de la raiz, el método entonces esta caracterizado por la sucesión:

\[ x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)} \]

Esto se generaliza a funciones \(f:\mathbb{R^n}\to\mathbb{R^m}\) (Método de Newton-Rhapson) con

\[ \vec{x_{i+1}}=\vec{x_i}-J(\vec{x_i})^{-1}f(\vec{x_i}) \]

Método de Newton

Ventajas

  • Es aplicable para cualquier función (difenciable con algunas restricciones en la segunda derivada)
  • Converge de manera cuadratica cuando se cumplen supuestos de suavidad de la función (Aproximadamente duplica los digitos exactos en cada iteración)

Desventajas

  • Puede no converger (converge en una vecindad de la solución).
  • Hay dificultad operacional al calcular la derivada (la matriz \(J\) y su inversa)
  • Encuentra solo una solución (en caso de existir varias)
  • No aprovecha la estructura del problema (optimización)

Paso 1: Dar un valor inicial a x.

Paso 2: Calcular la recta tangente de la función en el punto dado.

Paso 3: Usar la recta tangente para actualizar el x.

Paso 4: Repetir hasta la convergencia.

Visualización

Divergencia en el método de Newton

Método EM (Expectation Maximization)

En el ejemplo inicial tenemos una variable que no fue observada, de la cual solo podemos conjeturar su valor.

  • Es posible estimar la esperanza de esta variable si tuvieramos los valores de los parametros de media y varianza de las subpoblaciones.
  • Es posible estimar los parametros de las dos distribuciones normales si tuvieramos la probabilidad de pertenencia a cada grupo.

Pasos del método

  1. Inicializar el proceso con algun valor para los valores de los parametros a estimar.
  2. Optimizar la verosimilitud de los parametros ocultos considerando que los parametros visibles son conocidos (este es el paso E, ya que es equivalente a tomar la esperanza bajo los parametros desconocidos).
  3. Optimizar la verosimilitud de los parametros visibles considerando que los parametros ocultos son conocidos (Este es el paso M).
  4. Repetir hasta la convergencia

Representacion Gráfica

Drawing

Expectation Maximizaton

Ventajas

  • Siempre converge.
  • Relativamente intuitivo y facil de implementar.

Desventajas

  • Solo funciona para problemas de optimizacion que cumplen cierta estructura.
  • La velocidad de convergencia es baja.
  • Esta convergencia no se puede asegurar que sea a un optimo global.
  • Encuentra solo una solución (en caso de existir varias).

Para mixtura

En el caso especial de una mixtura (ejemplo inicial).

  • Inicie los parametros \(\mu_1,\sigma_1^2,\mu_2,\sigma_2^2.\pi\)
  • Expectation Step: calcule las probabilidates

\[ \gamma_i=\frac{\pi \Phi_{\mu_2}(y_i)}{\pi \Phi_{\mu_2}(y_i)+\pi \Phi_{\mu_1}(y_i)} \]

\[ \pi=\bar \gamma_i \]

  • Paso de Maximización: Obtenga el estimador maximo verosimil para \(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2\). Que son los estimadores usuales tomando como pesos las probabilidades \(\gamma_i\) calculadas en el paso anterior.

  • Repita los pasos dos y tres hasta la convergencia

Paso 0.

Se utilizaron números aleatorios entre 0 y 10 para inicializar las medias, las varianzas se iniciaron cada una con la mitad de la varianza muestral.

Paso 1.

Se actualiza la probabilidad de que cada punto pertenezca a cada poblacion (es decir se calcula \(\gamma_i\))

Paso 2.

Con esta información se estiman de nuevo medias y varianzas usando el estimador maximo verosimil.

Paso 3.

Repetir hasta la convergencia

Visualización

Referencias

  1. Hastie, T., Tibshirani, R.,, Friedman, J. (2001). The Elements of Statistical Learning. New York, NY, USA: Springer New York Inc.
  2. Kendall E. Atkinson, An Introduction to Numerical Analysis, (1989) John Wiley & Sons, Inc, ISBN 0-471-62489-6.
  3. Expectation Maximization, how it works. Consultado el 3 de mayo en https://www.youtube.com/watch?v=iQoXFmbXRJA.

Gracias

  • Esta presentación se realizó usando RMarkdown (ioslides)
  • Las gráficas de esta presentación se realizaron usando ggplot2
  • Las animaciones se realizaron usando plotly
  • Esta presentación y el código fuente están en https://github.com/kuhsibiris/Newton-y-EM bajo licencia creative commons